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Abstract 

In bacterial populations, cells are able to cooperate in order to yield complex collective func- 
tionalities. Interest in population-level cellular behaviour is increasing, due to both our expanding 
knowledge of the underlying biological principles, and the growing range of possible applications for 
engineered microbial consortia. Researchers in the field of synthetic biology - the application of en- 
gineering principles to living systems - have, for example, recently shown how useful decision-making 
circuits may be distributed across a bacterial population. The ability of cells to interact through 
small signalling molecules (a mechanism known as quorum sensing) is the basis for the majority of 
existing engineered systems. However, horizontal gene transfer (or conjugation) offers the possibility 
of cells exchanging messages (using DNA) that are much more information-rich. The potential of 
engineering this conjugation mechanism to suit specific goals will guide future developments in this 
area. Motivated by a lack of computational models for examining the specific dynamics of conju- 
gation, we present a simulation framework for its further study. We present an agent-based model 
for conjugation dynamics, with realistic handling of physical forces. Our framework combines the 
management of intercellular interactions together with simulation of intracellular genetic networks, 
to provide a general-purpose platform. We validate our simulations against existing experimental 
data, and then demonstrate how the emergent mixing patterns of multi-strain populations can affect 
conjugation dynamics. Our model of conjugation, based on a probability distribution, may be easily 
tuned to correspond to the behaviour of different cell types. Simulation code and movies are available 
at |http://code.google.com/p/discus/| 



1 Introduction 

Researchers in the emerging field of synthetic biol- 
ogy [T] have demonstrated the successful construc- 
tion of devices based on populations of engineered 
microbes [5J [3] 0] . Recent work has focussed atten- 
tion on the combination of single-cell intracellular 
devices [5J |5] with intercellular engineering, in order 
to build increasingly complex systems. 

To date, most work on engineered cell-cell com- 
munication has focussed on quorum-sensing (QS) 
[7], which may be thought of as a communication 
protocol to facilitate inter-bacterial communication 
via the generation and receiving of small signal 
molecules. However, recent studies on DNA mes- 
saging [8 highlight the importance and utility of 
transferring whole sets of DNA molecules from one 
cell (the so-called donor) to another (the recipient). 
Bacterial conjugation is a cell-to-cell communica- 
tion mechanism [HI [TO] that enables such transfers 
to occur. 

In this paper we present a simulation platform 
that realistically simulates (in a modular fashion) 



both intracellular genetic networks and intercellu- 
lar communication via conjugation. To our knowl- 
edge, this is the first such platform to offer both of 
these facilities. We first review previous work on 
cell simulation, before presenting the details of our 
model. We validate it against previous experimen- 
tal work, and then discuss possible applications of 
our method. 

2 Previous work 

The rapid development of bacterial-based devices 
is accompanied by a need for computational sim- 
ulations and mathematical modelling to facilitate 
the characterisation and design of such systems. A 
number of of platforms and methods are available 
for this purpose. Agent-based models (AbMs) are 
widely used [11], and were first used to study mi- 
crobial growth in BacSim |12j . Continuous models 
have also been proposed [13], and recent develop- 
ments make use of hardware optimisation, by using 
CPUs (Graphics Processing Units) in order to scale 



up the number of cells simulated [T3] . 

Because of the complexity of the system under 
study, several computational platforms focus on 
cither specific cellular behaviours (e.g., bacterial 
chemotaxis |15| . morphogenesis of dense tissue like 
systems |16j ) , or on specific organisms (e.g., Myx- 
ococcus xanthus |17j). Platforms that incorporate 
cell-cell communication generally focus their atten- 
tion on quorum-sensing. Simulations of conjugation 
do exist, but these consider cells as abstracted cir- 
cular objects [TS1 US] - We demonstrate in this paper 
how a consideration of the shape of cells is an es- 
sential feature for understanding the conjugation 
behaviour of the population. We now describe our 
model for bacterial growth, in which conjugation is 
handled explicitly. 



3 Methods 

We use an individual-based modelling approach [20] 
to the study of conjugation dynamics. This mod- 
els each cell as an individual, mobile entity, each 
of which is subject to physical forces arising from 
contact with other cells and the environment (e.g., 
surfaces). Each cell has a number of different at- 
tributes, listed in Table [T] which correspond to 
various physiological states and characteristics. 

Bacteria are modelled as rod-shaped cells with 
a constant radius (parameter width in Table [2]). 
Elongation processes occur along the longitudinal 
axis, which has a minimum dimension of length, 
and division takes place whenever a the cell mea- 
sures 2*length. The division of a cell into two new 
daughter cells is also controlled by max_overlap, 
which monitors the physical pressure affecting each 
cell; if the pressure exceeds this parameter value, 
the cell delays its growth and division. Thus, a cell 
with pressure grows slower than without it. In Fig- 
ure [T]^ we see an snapshot of a population with dif- 
ferent cell lengths, due to the pressure-dependent 
behaviour. The global parameter growth_speed 
(Table [2]) also helps us simulate cell flexibility in 
a realistic fashion. This parameter defines a "cut 
off" value for the number of iterations in which the 
physics engine must resolve all the current forces 
and collisions. Thus, smaller values will cause the 
solver to be effectively "overloaded" , and some colli- 



sions may, as a result, be partially undetected. This 
means that cells behave as flexible shapes, which 
gives the simulation a more realistic performance. 
In Figurejlp we show how changes in growth_speed 
affect the simulation, using bigger (left) to smaller 
(right) values. 

Horizontal genetic transfer (or conjugation) is 
modelled using an elastic spring to connect donor 
and recipient cells. Parameter c_time defines the 
duration of that linkage, which determines the time 
in which the DNA is transferred. The springs are 
constantly monitored to ensure that they physically 
connect both cells during conjugation. Importantly, 
during conjugation, the resolution of collisions in- 
volving relevant cells considers the forces produced 
by the spring connection, in order to calculate the 
final movement of the bacteria. By coupling cells in 
this way, we obtain realistic population-level phys- 
ical patterns that emerge as a result of large num- 
bers of conjugation events. Figure [Tj3 shows this 
process, with a donor cell (red) and a recipient cell 
(grey) which becomes a transconjugant (yellow). A 
transconjugant cell is one that was initially a re- 
cipient, but which has been conjugated during its 
lifetime. Thus, it already has the DNA information 
transferred by the donor. 

This agent-based algorithm has an iteration- 
driven structure, where - after initialisation of the 
main global parameters - it repeatedly performs the 
following steps for each cell: 



1. Update springs (position and timing). 

2. Perform cell division (if cell is ready). 

3. Elongate cell (every growth_speed steps). 

4. Handle conjugation. 

5. Update physical position. 

Conjugation decisions (step 4) made by cells are 
driven by three sequential steps: 

1. The cell decides, following a probability dis- 
tribution, whether or not to conjugate (one 
trial per iteration). 

2. If conjugating, randomly select a mate from 
surrounding bacteria (if present). 
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Table 1. Cell attributes. 



Attribute 


type 


Definition 


shape 


pymunk.Sheipc 


Shape of the cell 


program 


[m ... mj 


List of the i regulatory network molecules (m) 


elongation 


[int, int] 


Elongation values (one per cell pole) 


position 


[x,y] 


Coordinates of centre point, x and y 


speed 


float 


Velocity 


conjugating 


Boolean 


Conjugation state 


plasmid 


Boolean 


Program state (present /not present) 


role 


int 


Donor (0), recipient (1) or transconjugant (2) 


partner 


int 


Role of plasmid transfer cell 




Overlapping + 



Figure 1. Cell behaviours at low scale. (A): Different cell length due to asynchronous growth (pressure 
dependent). Two cells marked with red arrows. (B): A donor cell (red) starts the conjugation process 
with a recipient (grey) which turns into transconjugant (yellow). The pilus (green) is an elastic spring 
that links the two cells until the process is finished. (C): Different overlapping levels within the cells of 
a population. 



3. If valid mate is found, effect conjugation 
transfer. 

The discrete probability distribution used for 
the conjugation process is C(N, p, c_time), where 
N is the number of trials in a cell lifetime (width 
* length) , p is the success probability in each trial 
(with p G [0. . .1]) and c_time is the time interval 
during p = 0.0 (i.e., when the cell is already con- 
jugating). As stated in [21], p can vary, depending 
on whether the cell is a donor (p_d), a transconju- 
gant that received the DNA message from a donor 
(p_tl), or a transconjugant that received the DNA 
from another transconjugant (p_t2). 

Intracellular circuitry is modelled separately, 
and then introduced into each cell by storing the 
state of the circuit in an attribute of the cell 
(program) . Thus, there are effectively as many 
copies of the circuit as cells in the simulation. This 
circuit simulation is implemented in a modular fash- 
ion, so that the internal cellular "program" may be 
easily replaced with any other. In this paper we 
demonstrate the principle using a two-component 
genetic oscillator as the DNA message that is ex- 
changed through conjugation. The ordinary differ- 
ential equations (ODEs) for this circuit are: 

dx ( 1 + ax 2 \ 
dt \ 1 + x A + oy J 

dy 1 + ax 2 

dt 1 + x A 

which are detailed in [35], as well as the meaning 
and value of the parameters (we use the same values 
in the code provided). We have also recently used 
our software platform to investigate the spatial be- 
haviour of a reconfigurable genetic logic circuit |23| , 
which demonstrates how it may easily be modified 
to accommodate different sets of equations. The ac- 
tions controlling the growth rates of cells occur on 
a longer time scale than the integration steps that 
control molecular reactions (as equations [I] and [2]) . 
In order to ensure synchronisation, the parameter 
network_steps defines the number of integration 
steps of the ODEs that run per Gt. Thus, a number 
of network_steps/Gt integration steps will update 
the attribute network of each cell every iteration. 

Other important physical parameters listed in 
Table [Hare 

spring_rest_length, spring_stif f ness and 



spring_damping; these are three parameters to" 
model the material and behaviour of the bacterial 
pilus (i.e. the spring) during conjugation. Pa- 
rameter cell_infancy is a delay period, during 
which a cell is considered to be too young to conju- 
gate (as observed experimentally [2T] ). Parameters 
pymunk_steps and pymunk_clock_ticks are used 
by the physics engine to update the world, and may 
be adjusted by the user in order to alter the per- 
formance of the simulation (machine dependent). 
Parameters bacjnass and bac_f riction play a role 
in collision handling. 

Our platform is writting in Python, and 
makes use of the physics engine pymunk 
(www.pymunk.org) as a wrapper for the 2D 
physics library Chipmunk, which is written in C 
(www.chipmunk-physics.net/). As cells are rep- 
resented as semi-rigid bodies in a 2D lattice, py- 
munk handles the physical environment on our be- 
half. For monitoring purposes, parameters Gt and 
real_Gt allow us to stablise the relation between it- 
erations and clock minutes: minute = Gt/real_Gt 
(units: iterations). The platform, which we call 
DiSCUS (Discrete Simulation of Conjugation Us- 
ing Springs) is available at the project repository 
at |http://code. google. com/p / discus / 



4 Results 

We now describe the results of experiments to vali- 
date our conjugation model, using four sets of sim- 
ulations. As we aim to understand the behaviour 
of cells in small-scale two-dimensional populations 
(as occur in microfluidic environments), we avoid 
the sorts of extreme overlapping situations shown 
in Figure [lp (right). We first validate individual 
conjugation dynamics; then we validate the biome- 
chanical properties of the simulation; the third set 
of experiments concerns the transfer of the two- 
component oscillator, and the final set of experi- 
ments study the effects of mixing on conjugation 
dynamics. We now describe in detail the results of 
each set of experiments. 

4.1 Conjugation dynamics 

The objective of the first set of experiments is to val- 
idate the software in terms of conjugation dynam- 
ics. For that purpose, we first focus on conjugation, 
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Figure 2. Validation of cell movement and conjugation dynamics using real data. (A): Figure 
extracted from [5T] where a colony of Pseudomonas putida is divided into dark red donor cells (DsRed), 
yellow recipient cells (YFP) and transconjugants, expressing both yellow and green light (YFP and 
GFP). The upper row shows the transconjugant signal, and the bottom row shows the whole 
community. (B and C): Simulation results. Two simulations of similar colonies are recorded over 
exactly the same time intervals (min). The colours of the cells match the colours observed in (A). 
Graphs (D), (F) and( H) are extracted from |24j . and show experimental results of Escherichia Coli 
growth regarding density, velocity and ordering (respectively). Graphs (E), (G) and (I) correspond to 
simulations in similar conditions to [23] , for the same parameters (density, velocity and ordering 
respectively). Tests 1, 2 and 3 in graphs correspond to different spatial distribution of cells inside the 
microfluidic chanel (details in text). 
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Figure 3. Horizontal transfer of a two-component genetic oscillator. A: Transcriptional-level design. 
The activator (green) acts on itself and on the repressor (red) by inducing the transcription of both. 
The represor acts on the activator by repressing its transcription. As a result, molecule x (as well as y 
oscillates in time. B: Cells growing in a cross-shaped channel. At 250 (minutes) we clearly see the 
position of both donor (left hand, with the oscillator inside) and recipient (right hand, empty) strains. 
The intensity of green colour denotes the amount of molecule x inside cells in that specific time. 
Through conjugation, the oscillator is copied to the initial recipient strain (t = 1000 min). C: Using the 
same experiment as in B, measurement over time of the maximum level of molecule a; in a single cell. 
Black arrow highlights the point in time when conjugation starts (t ~ 550 min). 



using images of a Pseudomonas putida population 
(Figure [2|\_) extracted with permission from [2"T] , 
These show donor cells (dark red) growing in con- 
tact with recipients (yellow). The DNA information 
they share after conjugation makes the transconju- 
gant cells display GFP (green fluorescent protein). 
We adjusted the parameters of our simulations until 
the behaviour matched the images of real cells (two 
simulations shown: Figures[2j3 and[2jll), in terms of 
both time-series behaviour and the type of physical 
pattern displayed. It is important to note that the 
differential probabilities of conjugation of donors 
and transconjugants (higher in the latter) causes di- 
rectional spreading of the DNA information. After 
the first transconjugant appears (160 minutes), the 
newly-formed transconjugants appear -most prob- 
ably - in the immediate neighbourhood. The final 
parameter values used to reproduce this experiment 
are: width=5, length=15, growth_speed=30, 
p_d=0.001, p_tl=0.02, p_t2=0.05 and c_time=450 
(the rest of the parameters are as defined in the 
DiSCUS distribution). Movie DemoConjugationl 
(found in the project repository) shows a simulation 
of a similar experiment where the transconjugants 
do not act as new donors. 

4.2 Biomechanical properties 

The second set of validation experiments focuses on 
biomechanical movement. We use data from |24j . 
which describe an Escherichia coli colony growing 
in a microfluidic channel (30 * 50 * 1 /mi 3 ) (Fig- 
ures [2p, [2^ and [2jT) . Using the same parameter 
setup (width=5, length=24, growth_speed=30) 
we highlight how different initial positioning of cells 
inside the channel can affect the final result (testl, 
with more cells observed in the centre than at the 
edges; testB with all cells initially in the centre; testS 
with all cells homogeneously spread along the chan- 
nel). Density graphs (Figures [2p and[2j3) show the 
increasing curve as the channel becomes more pop- 
ulated (results vary depending on which area is con- 
sidered for monitoring). Velocity gradients (Figures 
[2^ and [2]GI) depict the differential velocity across 
the longitudinal axis of the channel with respect to 
the centre (we see negative values when the cells in 
the centre move faster than the rest). The differ- 
ence in the y axis is due to our considering different 
spacial intervals in the velocity gradient calculation. 
Ordering graphs (Figures [2ji and[2j) are based on 



calculating the cosine of a cell's angle with respect 
to the longitudinal axis of the channel (e.g. angle 
0, cos(0)=l, completely aligned). As time increases, 
we see that the cells tend to align themselves. 

4.3 Internal cell "program" 

Figure [3] shows the results of the next set of exper- 
iments, aiming at studying the horizontal transfer 
of a two-component oscillator. The transcriptional- 
level design of the circuit is shown in Figure ^ A) , 
which causes the molecular concentration of a re- 
pressor (y) and an activator (x) to oscillate over 
time. Each molecule is produced by a gene when 
its upstream promoter is activated. The activator 
can induce its own production at the same time as 
inducing the production of the repressor, which in 
turns inhibits the production of the activator. In 
Figure [3^3 we place (250 minutes) an initial donor 
colony (with the two-component oscillator inside) 
on the left-hand side of a cross-shaped channel, 
while a recipient colony is placed on the right. 

As the equations for the oscillator have no 
stochasticity, every cell of the donor strain shows 
exactly the same state of the circuit as every other 
cell. At the beginning (250 minutes), these cells 
show a green colour (corresponding to the molecular 
concentration of x) which is switched off during the 
time intervals in which the repressor is on (see time 
profile of Figure j3p). When conjugation starts (at 
around t ~ 550), the newly formed transconjugant 
cells are given the circuit but, importantly, they do 
not share the state of the circuit of the cells from 
which they receive the message. During the DNA 
transfer, it is only the plasmid (circuit carrier) that 
is copied into the recipient; therefore, both molec- 
ular concentrations are null, and the circuit begins 
its functioning from the initial stage. That is why 
we clearly observe different green intensities within 
the community. This asynchronous behaviour hap- 
pens only in the transconjugants, while the circuits 
inside the donors always run synchronously (due to 
deterministic equations) . 

A time profile of the previous experiment is 
shown in Figure [3p, where the maximum level of 
activator (x) concentration in a single cell (com- 
pared with the whole population) is recorded over 
time. Before conjugation starts, all cells in the con- 
sortia display perfect synchrony. After conjugation 
(shown with an arrow on the graph) there is al- 



ways a cell with the maximum level of activator, 
which demonstrates high asynchrony. All param- 
eter values regarding cell dimensions or conjuga- 
tion probabilities are the same as in Figures [2j3 
and [2fl Parameters relevant to the oscillator are: 
network_steps=18 and Gt=450. Movie DemoCon- 
jutagion2 (found in the project repository) shows 
this experiment and Movie D emo Dynamics 1 shows 
donor cells growing with a stochastic version of the 
oscillator. 

4.4 Effects of mixing 

Conjugation behaviour within a population may be 
altered in different ways to achieve different be- 
haviours, depending on the desired application. For 
example, in the previous experiments described in 
this paper, transconjugants are unable to act as re- 
cipients (simulating a radical entry exclusion |25|). 
That is to say, they will not receive more plasmids 
(genetic circuits) from either donors or transcon- 
jugants. Furthermore, we may also engineer the 
transconjugants to stop acting as new donors [TU] . 
so that only the original donors have the ability 
to transfer the DNA message. Mixing of the cell 
population becomes essential in this last scenario, 
in order to ensure maximal contact between donors 
and recipients. In Figure [3j3 we see how, at the end 
of the experiment (fOOO minutes), donors and the 
transconjugants cover different areas of the channel 
(left and right respectively), without being mixed. 

We now study the autonomous mixing be- 
haviour of cells under different environmental condi- 
tions, with the third sets of experiments (Figure [4]). 
Firstly we investigate how morphological changes 
in a longitudinal microfluidic channel can affect the 
patterns being formed by the consortia and its mix- 
ing. Figure |4]A. shows three bacterial strains (each 
shown in a different colour) growing in a channel 
from different starting points. As we can see, their 
mixing is highly improbable. The main reason for 
this is the velocity and directionality of the cells. As 
the cells are washed out at the edges of the channel, 
all of them travel (they only have passive movement 
while being pushed) at variable speed from centre 
to left or from centre to right. This causes the cells 
to have the same direction (see vector field), which 
in turn makes mixing more difficult . In Figures [4j3 
and [4p we show the result of altering the morphol- 
ogy of the channels, by adding columns and zig-zag 



walls, respectively. As a result, the cells show differ? 
cnt directionality (see vector fields) and the strains 
have a higher probability of being mixed. In both 
experiments (unlike Q A)), the three strains are in 
contact at some point. If the experimental applica- 
tion relied on the conjugation of purple and yellow 
cell pairs, for example, we can see that these phys- 
ical changes in the channel would be essential. 

Another way to intensify the mixing of strains 
in a microfluidic trap is to change the main channel 
flow strength, with the objective of creating turbu- 
lence inside the trap. Figures [4jZ) and |4j5 show a 
three-strain population growing in a trap (identi- 
cal initial positioning of cells in both) with the only 
difference being that the strong flow in the main 
channel (white arrows) of creates turbulence (in 
the direction of the circled arrow). Two different 
colonies are simulated in each experiment, inside 
both symmetrical and independent traps. We see 
how turbulence helps the cells to get mixed, thanks 
to the constant change of direction they display 
(see vector field in Figure]^). Furthermore, we can 
avoid missing one strain, as happens with the yel- 
low cells in [4jZ> (see Movie DemoDynamics2 -found 
in the project repository-), where another run of 
this experiment is shown). 

Investigations of how manual mixing can affect 
conjugation frequencies are described in in [10], us- 
ing an Escherichia coli population. We now repro- 
duce those results using our software, and give valu- 
able insight into the reasons for that behaviour: the 
isolation of the recipients. For that purpose (Figure 
[5]) we grow a population of donors (D, red) and re- 
cipients (R, yellow) in which the ratio D/R is 50% 
and the transconjungants (T, green) are unable to 
act as new donors. The frequency of conjugation, 
Y, is measured as Y = T/(R + T). The graph in 
Figure [5p shows the frequency after 560 minutes of 
untouched populations (not mixed, blue bars) and 
populations that have been manually mixed at 420 
minutes (red bars). The difference that the mix- 
ing produces is based on the isolation of the recip- 
ients in untouched populations. Figure [5]A. shows 
two different occasions in which clusters of recipi- 
ents are formed, where the transconjugants do not 
allow donors to reach new possible mates. After the 
population is completely "shuffled" ([5^3), the clus- 
ters are dissolved, and new pairs of donor-recipient 
can arise in the new topology. 





Figure 4. Dynamical mixing of bacterial strains under different conditions. (A): Three strains (purple, 
yellow and blue) growing in a longitudinal pipe. Detail (on the right) shows the vector field 
corresponding to the red square (on the left) where the arrows display both directionality and velocity 
of every cell. (B): Similar to (A), but with four columns placed in the middle of the pipe. (C): Similar 
to A, but with zig-zag borders along the pipe. In (A), (B) and (C), the speed in vector fields is 
measured in arbitrary units (a.u.). (D): Six strains (three per trap) grow in two square traps on one 
side of the main longitudinal channel. The flow in the channels follows the direction of the arrows. (E): 
Similar to (D) but with a much stronger flow in channels, causing turbulence in traps (long circled 
arrow). (F): Vector field of experiment (E). 
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Cell dimensions 

Figure 5. Effects of manual mixing on conjugation frequency. (A): Recipient-trapping behaviour of a 
population with donors (red), transconjugants (green) and recipients (yellow). Two snapshots depict 
clearly-observed clusters. (B): Population after random mixing, where the clusters are automatically 
dissolved. (C): Graph showing conjugation frequencies (Y = T/(R + T)) of 560-minute experiments 
(ratio D/R = 50%). Blue bars represent Y on an untouched population, while red bars represent Y 
when the population is mixed at 420 minutes. The two sets of bars correspond to experiments with 
different cell dimensions (1x3 -left- and 1x2 -right). Error bars show variation across 15 experiments of 
each class. 




An interesting result from Figure [5p is the fact 
that the smaller the size of the cell, the higher re- 
sults we observe for conjugation frequencies. This 
may be due to the fact that smaller cells are able to 
slip through physical gaps, and the biomechanical 
ordering of the population becomes more "fuzzy" . 
This underlines the importance of considering the 
physical shape of cells, since circle-shaped cells 
would not give valid results. 

5 Discussion 

The conjugation model presented here is the first 
agent-based model to explicitly simulates the con- 
jugation process with growing rod-shaped cells. Full 
validation against real data is performed, which 
shows the capacity of the software to reproduce ob- 
served behaviour. In addition, the mixing study of- 
fers valuable insights into the design of multi-strain 
populations. The software also allows for genetic 
programs to be installed inside cells; the potential 
for horizontal gene transfer to recreate distributed 
information processing within a microbial consor- 
tium is of significant interest in synthetic biology, 
and the software presented will aid the design and 
testing of systems before their wet-lab implementa- 
tion. 
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